rm(list = ls())
library(foreign)
library(fwb)
library(dplyr)
library(brglm2)
library(beepr)
library(modelsummary)
library(lmtest)
library(sandwich)
library(marginaleffects)
library(kableExtra)
library(marginaleffects)
library(ggplot2)
library(ggbreak)

options(modelsummary_output = "latex_tabularray")
# Remove the current version first (optional, but recommended)
remove.packages("modelsummary")
set.seed(555)

# Install the older version from the CRAN archives
remotes::install_version("modelsummary", version = "0.10.0")
shame2 = read.csv("/Users/cbrandt/Dropbox/AAA_Projects/UNSC/R code/replication_data/brandt_kotajoki_naming_and_shaming_at_the_UNSC_panel_data.csv")

universal_r = 1000
shame2$e10 = ifelse(shame2$p5==1, 0, 1 )
shame2$shamed_dummy <- as.integer(shame2$shamed_dummy)

#Appendix Summary statistics table -----
shame2_clean <- shame2[, c("shamed_dummy", "p5", "v2x_clphy_lag", "best_fatality_estimate_lag_log", 
                           "ext_confirmed_lag", "battledeaths_lag_log", 
                           "incompatibility", "pk_lag_agg_log", "any_member_shamed_year_prior", 
                           "total_resolutions_log_lag")]


shame2_clean <- shame2_clean[complete.cases(shame2_clean), ]
stargazer(as.data.frame(shame2_clean), title="Descriptive Statistics",
          covariate.labels=c("Naming and Shaming", "P5", "Physical Integrity",
                             "One-Sided Violence","External Support", "Intensity",
                             "Secessionist","Peacekeeping", "Prior Shame", "Resolutions", "Islam", "Post 9/11"),
          out="/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/appendix_summary_statistics.tex")


cm  = c("cw"=  "Civil War",  "isfighting3"="Fighting NSA",
        'ext_confirmed_lag' = 'External Support',
        'ext_confirmed_lag:p5' = 'External Support*P5',
        'p5:ext_confirmed_lag' = 'External Support*P5',
        'best_fatality_estimate_lag_log' = 'One-Sided Violence',
        'best_fatality_estimate_lag_log:p5' = 'One-Sided Violence*P5',
        'p5:best_fatality_estimate_lag_log' = 'One-Sided Violence*P5',
        'ext_confirmed_lag:best_fatality_estimate_lag_log' = 'Ext. Sup.*OSV',
        'best_fatality_estimate_lag_log:ext_confirmed_lag' = 'Ext. Sup.*OSV',
        "no_of_allies_fighting_rebel_dummy"= "Allies",
        "no_of_allies_fighting_rebel_dummy:p5"= "Allies*P5",
        "p5:no_of_allies_fighting_rebel_dummy"= "Allies*P5",
        'ext_confirmed_lag:best_fatality_estimate_lag_log:p5' = 'Ext. Sup.*OSV*P5',
        'best_fatality_estimate_lag_log:ext_confirmed_lag:p5' = 'Ext. Sup.*OSV*P5',
        "best_fatality_estimate_lag_log:no_of_allies_fighting_rebel_dummy" = "Allies*OSV",
        'no_of_allies_fighting_rebel_dummy:best_fatality_estimate_lag_log' = 'Allies*OSV',
        
        "ext_nsa_confirmed_lag"="Ext Support to NSAG",
        "ext_nsa_confirmed_lag:p5"="Ext Support to NSAG*P5",
        "p5:ext_nsa_confirmed_lag"="Ext Support to NSAG*P5",
        
        "post911"="post911",
        "any_osv" = "Any OSV",
        "p5" = "P5", 
        "v2x_clphy_lag" = "Physical Integrity",
        "v2x_clphy_lag:best_fatality_estimate" = "One-Sided Vio. x Phys. Int.",
        "any_osv:ext_confirmed" = "Ext Support*Any OSV",
        
        "polity2" = "Regime Type",
        "battledeaths_lag_log"= "Intensity",
        "strength" = "Strength",
        "incompatibility"= "Secessionist",
        "hrw_prev_2" = "SV: HRW2",
        "state_prev_2" = "SV: State2",
        "ai_prev_2" = "SV: AI2",
        "rebpolwing" = "Political Wing",
        "pk_lag_agg_log" = "Peacekeeping",
        
        
        "prior_no_russia" = "Prior Shame*",
        "prior_no_china" = "Prior Shame*",
        "prior_no_us" = "Prior Shame*",
        "prior_no_uk" = "Prior Shame*",
        "prior_no_france" = "Prior Shame*",
        "prior_no_russia_china" = "Prior Shame*",
        "prior_no_western" = "Prior Shame*",
        
        
        "any_member_shamed_year_prior"= "Prior Shame",
        "total_resolutions_log_lag" = "Resolutions",
        "t_a"="Terrorism1",
        "t_b"="Terrorism2",
        "k_a"="Terrorism3",
        "k_b"="Terrorism4",
        
        "us_shamed_same_year_prior"= "US Shamed Prior",
        "russia:us_shamed_same_year_prior" = "US Shamed Prior*Russia Fe",
        "us_shamed_same_year_prior:russia" = "US Shamed Prior*Russia Fe",
        "china:us_shamed_same_year_prior" = "US Shamed Prior*China FE",
        "us_shamed_same_year_prior:china" = "US Shamed Prior*China FE",
        "uk_shamed_same_year_prior"= "UK Shamed Prior",
        "france_shamed_same_year_prior"= "France Shamed Prior",
        "russia_shamed_same_year_prior"= "Russia Shamed Prior",
        "china_shamed_same_year_prior"= "China Shamed Prior",
        
        "western" = "Western",
        "russ_china" = "Rus/China",
        "western_shame_prior" = "Western Shamed Prior", 
        "russ_china_shame_prior" = "Russia/China Shamed Prior",
        
        "western_shame_prior:russ_china" =  "Western Prior*RC", 
        "western_shame_prior:russia" =  "Western Prior*Russia", 
        "western_shame_prior:china" =  "Western Prior*China", 
        "western_shame_prior:western" =  "Western Prior*Western",
        "russ_china_shame_prior:western" = "RC Prior*Western",
        "p5_shamed_same_year_prior"= "P5 Shamed Prior",
        "p5:p5_shamed_same_year_prior" = "P5 Shamed Prior*P5",
        "e10_shamed_same_year_prior"= "E10 Shamed Prior",
        "e10_shamed_same_year_prior:p5" = "E10 Shamed Prior*P5",
        "rolling_shaming_group_sum_5yr" = "Rolling by Group", 
        "rolling_shaming_state_sum_5yr" = "Rolling by State",
        "nego" = "Negotiations",
        "usa" = "US FE",
        "russ_china_shame_prior:usa" = "Russia/China Prior*US FE",
        
        "uk" = "UK FE",
        "russ_china_shame_prior:uk" = "Russia/China Prior*UK FE",
        "france" = "France FE",
        "russ_china_shame_prior:france" = "Russia/China Prior*Fr FE",
        "russia" = "Russia FE",
        "china" = "China FE",
        
        "usa:russ_china_shame_prior" = "Russia/China Prior*US FE",
        "uk:russ_china_shame_prior" = "Russia/China Prior*UK FE",
        "france:russ_china_shame_prior" = "Russia/China Prior*Fr FE",
        
        "us_shamed_same_year_prior:usa" = "US Shamed Prior*US FE",
        "usa:uk_shamed_same_year_prior" = "UK Prior*US FE",
        "usa:france_shamed_same_year_prior" = "France Prior*US FE",
        "usa:russia_shamed_same_year_prior" = "Russia Prior*US FE",
        "usa:china_shamed_same_year_prior" = "China Prior*US FE",
        
        "uk:us_shamed_same_year_prior" = "US Shamed Prior*UK FE",
        "uk:uk_shamed_same_year_prior" = "UK Prior*UK FE",
        "uk:france_shamed_same_year_prior" = "France Prior*UK FE",
        "uk:russia_shamed_same_year_prior" = "Russia Prior*UK FE",
        "uk:china_shamed_same_year_prior" = "China Prior*UK FE",
        
        "france:us_shamed_same_year_prior" = "US Shamed Prior*Fr FE",
        "france:uk_shamed_same_year_prior" = "UK Prior*Fr FE",
        "france:france_shamed_same_year_prior" = "France Prior*Fr FE",
        "france:russia_shamed_same_year_prior" = "Russia Prior*Fr FE",
        "france:china_shamed_same_year_prior" = "China Prior*Fr FE",
        
        
        "russia:uk_shamed_same_year_prior" = "UK Prior*Ru FE",
        "russia:france_shamed_same_year_prior" = "France Prior*Ru FE",
        "russia:russia_shamed_same_year_prior" = "Russia Prior*Ru FE",
        "russia:china_shamed_same_year_prior" = "China Prior*Ru FE",
        
        "china:uk_shamed_same_year_prior" = "UK Prior*Ch FE",
        "china:france_shamed_same_year_prior" = "France Prior*Ch FE",
        "china:russia_shamed_same_year_prior" = "Russia Prior*Ch FE",
        "china:china_shamed_same_year_prior" = "China Prior*Ch FE",
        "ally_shamed_prior_2_years" = "ally Shamed Prior",
        'strength_lag' = "Strength",
        "t_a_lag"="Terrorism1",
        "t_b_lag"="Terrorism2",
        "k_a_lag"="Terrorism3",
        "k_b_lag"="Terrorism4",
        'neg_lag' = 'Negotiations',
        
        'islam' = "Islamist",
        'post_911' = "Post 9/11",
        'islam:post_911' = "Islamist*Post 9/11",
        
        "year" = "Time",
        
        "e10_shamed_same_year_prior:best_fatality_estimate_lag_log" = "E10 Shamed Prior*OSV",
        
        "best_fatality_estimate_lag_log:p5_shamed_same_year_prior" = "P5 Shamed Prior*OSV",
        
        "e10_shamed_same_year_prior:best_fatality_estimate_lag_log:p5" = "E10 Shamed Prior*OSV*P5",
        
        "best_fatality_estimate_lag_log:p5:p5_shamed_same_year_prior" ="P5 Shamed Prior*OSV*P5", 
        
        "e10_shamed_same_year_prior:ext_confirmed_lag" = "E10 Shamed Prior*Ext. Sup.",
        
        "ext_confirmed_lag:p5_shamed_same_year_prior" = "P5 Shamed Prior*Ext. Sup.",
        "e10_shamed_same_year_prior:ext_confirmed_lag:p5" = "E10 Shamed Prior*Ext. Sup.*P5",
        "ext_confirmed_lag:p5:p5_shamed_same_year_prior" = "P5 Shamed Prior*Ext. Sup.*P5",
        "(Intercept)"="Intercept"
)

#####Appendix Tables =======
#Rare events logit-----
formula1 <- shamed_dummy ~ ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior

formula2 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*p5  +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag

formula3 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag - 1

# The argument type = "AS_mean" specifies Firth's method.
# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))

data_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

model1_br <- glm(shamed_dummy ~ ext_confirmed_lag * p5 +battledeaths_lag_log + v2x_clphy_lag + 
                   incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior, 
                 data = data_clean_all, 
                 family = binomial(link = "logit"),
                 method = brglm2::brglmFit,
                 type = "AS_mean")

model2_br <- glm(shamed_dummy ~ best_fatality_estimate_lag_log*p5+battledeaths_lag_log + v2x_clphy_lag + 
                   incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior, 
                 data = data_clean_all, 
                 family = binomial(link = "logit"),
                 method = brglm2::brglmFit,
                 type = "AS_mean")

model3_br <- glm(shamed_dummy ~ any_member_shamed_year_prior + 
                   best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
                   incompatibility + pk_lag_agg_log + total_resolutions_log_lag, 
                 data = data_clean_all, 
                 family = binomial(link = "logit"),
                 method = brglm2::brglmFit,
                 type = "AS_mean")

clustered_se1 <- coeftest(model1_br, vcov = vcovCL, 
                          cluster = ~unsc_gwn, data = data_clean_all)[,2]

clustered_se2 <- coeftest(model2_br, vcov = vcovCL, 
                          cluster = ~unsc_gwn, data = data_clean_all)[,2]
clustered_se3 <- coeftest(model3_br, vcov = vcovCL, 
                          cluster = ~unsc_gwn, data = data_clean_all)[,2]

calculate_mcfadden <- function(full_model) {
  # Get the log-likelihood of the full model
  loglik_full <- logLik(full_model)
  
  # Manually create the call for the null model. This is more robust
  # than update() for special model types like 'brglmFit'.
  null_call <- getCall(full_model)
  
  # Replace the original formula in the call with the null formula (e.g., "shamed_dummy ~ 1")
  # reformulate() is a safe way to build formulas
  null_call$formula <- reformulate("1", response = all.vars(formula(full_model))[1])
  
  # Now, execute this new call to fit the null model correctly
  # The 'envir' argument helps R find the data used in the original model
  null_model <- eval(null_call, envir = parent.frame())
  
  loglik_null <- logLik(null_model)
  
  # Calculate McFadden's R-squared
  mcfadden_r2 <- 1 - (loglik_full / loglik_null)
  
  return(as.numeric(mcfadden_r2))
}

mcfadden_rsquared1 = calculate_mcfadden(model1_br)
mcfadden_rsquared2 = calculate_mcfadden(model2_br)
mcfadden_rsquared3 = calculate_mcfadden(model3_br)

row <- data.frame("Coefficients" = "McFadden’s R2",
                  "Model 1" = mcfadden_rsquared1,
                  "Model 2" = mcfadden_rsquared2,
                  "Model 3" = mcfadden_rsquared3)
attr(row, "position") <- 26


msummary(
  list("Model 1" = model1_br, "Model 2" = model2_br, "Model 3" = model3_br),
  vcov = list(clustered_se1, clustered_se2,clustered_se3), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c("*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  # Customize the bottom section of the table
  
  #gof_custom = custom_gof,
  add_rows=row,
  title = "Odds Ratios from Rare-Events Logistic Regression on \\textit{Naming and Shaming}", 
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/rare_events_logit.tex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')
)

#IV = formal alliance -----
formula1 <- shamed_dummy ~ no_of_allies_fighting_rebel_dummy*p5 +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior
formula2 <- shamed_dummy ~ best_fatality_estimate_lag_log*p5 +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior
formula3 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*p5 + no_of_allies_fighting_rebel_dummy*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag

# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn", "p5"))

data_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()


# This code should now run without error
model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))

vcov_matrix1 <- vcovFWB(model1, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = data_clean_all$unsc_gwn, R = universal_r)

msummary(
  # The nested list for models is CORRECT for creating headers
  list(model1, model2, model3), 
  #models_with_headers, # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2,vcov_matrix3),
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c( "*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE,  
  title =  "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming}: Operationalizes Partnerships as Formal Alliances", 
  output = "latex",
  
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')) %>%
  #add_header_above(header = c(" " = 1, "E10" = 3, "P5" = 3)) %>%
  save_kable(file = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/appendix_formal_alliance.tex")

#Interaction between OSV and External Support -----
formula1 <- shamed_dummy ~ ext_confirmed_lag*best_fatality_estimate_lag_log+battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior
formula2 <- shamed_dummy ~ ext_confirmed_lag*best_fatality_estimate_lag_log*p5+battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior

# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))

data_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

# This code should now run without error
model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))


vcov_matrix1 <- vcovFWB(model1, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = data_clean_all$unsc_gwn, R = universal_r)


msummary(
  list(model1, model2), # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c("*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  # Customize the bottom section of the table
  title =  "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming}: Includes Interaction Between Strategic Partnerships and Protection Norms", 
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/appendix_interaction_osv_ext_suprt.tex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')
)

#Controls for rebel group strength-------
formula1 <- shamed_dummy ~ ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior
formula2 <- shamed_dummy ~ best_fatality_estimate_lag_log*p5 +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior
formula3 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + strength_lag


# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))

data_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

# This code should now run without error
model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))

vcov_matrix1 <- vcovFWB(model1, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = data_clean_all$unsc_gwn, R = universal_r)


# Generate the table correctly
msummary(
  list(model1, model2, model3), # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2,vcov_matrix3), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c("+" = .1, "*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  # Customize the bottom section of the table
  title = "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming}: Includes Rebel Group Strength", 
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/includes_rebel_group_strength.tex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')
)

# Interdependence between naming and shaming ----
formula1 <- shamed_dummy ~   
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + e10_shamed_same_year_prior
formula2 <- shamed_dummy ~   
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + p5_shamed_same_year_prior
formula3 <- shamed_dummy ~   
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + p5_shamed_same_year_prior*e10_shamed_same_year_prior

# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))

data_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

# This code should now run without error
model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))

vcov_matrix1 <- vcovFWB(model1, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = data_clean_all$unsc_gwn, R = universal_r)


# Generate the table correctly
msummary(
  list(model1, model2, model3), # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2,vcov_matrix3), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c("*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  title= "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming}: Includes P5 and E10 Prior Shaming",
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/interdependence1.tex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')
)

#Interdependence - separating out each P5 state -----
formula1 <- shamed_dummy ~ 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ us_shamed_same_year_prior + prior_no_us

formula2 <- shamed_dummy ~ 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ uk_shamed_same_year_prior + prior_no_uk

formula3 <- shamed_dummy ~ 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + france_shamed_same_year_prior + prior_no_france

formula4 <- shamed_dummy ~ 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + russia_shamed_same_year_prior + prior_no_russia

formula5 <- shamed_dummy ~ 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ china_shamed_same_year_prior + prior_no_china

formula6 <- shamed_dummy ~   
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + us_shamed_same_year_prior + 
  uk_shamed_same_year_prior + france_shamed_same_year_prior+ russia_shamed_same_year_prior+
  china_shamed_same_year_prior + e10_shamed_same_year_prior


# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3, formula4, formula5, formula6)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))

data_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))
model4 <- glm(formula4, data = data_clean_all, family = binomial("logit"))
model5 <- glm(formula5, data = data_clean_all, family = binomial("logit"))
model6 <- glm(formula6, data = data_clean_all, family = binomial("logit"))

vcov_matrix1 <- vcovFWB(model1, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix4 <- vcovFWB(model4, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix5 <- vcovFWB(model5, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix6 <- vcovFWB(model6, cluster = data_clean_all$unsc_gwn, R = universal_r)

# Generate the table correctly
msummary(
  list(model1, model2, model3, model4, model5, model6), # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2,vcov_matrix3, vcov_matrix4, vcov_matrix5,vcov_matrix6), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c("*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  title= "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming}: Includes Individual P5 Prior Shaming",
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/interdependence2.tex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')
)

#Interdependence -Russia/China v Western-----
model1 = model2 = model3 = model4 =NULL 
vcov_matrix1 = vcov_matrix2 = vcov_matrix3 = vcov_matrix4 = NULL 

formula1 <- shamed_dummy ~ 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+  us_shamed_same_year_prior*russia + prior_no_us

formula2 <- shamed_dummy ~ 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ us_shamed_same_year_prior*china + prior_no_us

formula3 <- shamed_dummy ~ 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + western_shame_prior*russia + prior_no_western

formula4 <- shamed_dummy ~ 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + western_shame_prior*china + prior_no_western

# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3, formula4)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))

data_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))
model4 <- glm(formula4, data = data_clean_all, family = binomial("logit"))


vcov_matrix1 <- vcovFWB(model1, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix4 <- vcovFWB(model4, cluster = data_clean_all$unsc_gwn, R = universal_r)


# Generate the table correctly
msummary(
  list(model1, model2, model3, model4), # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2,vcov_matrix3, vcov_matrix4), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c("*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  # Customize the bottom section of the table
  title= "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming}: Do Russia and China Respond to Western Powers?",
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/interdependence3.tex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')
)


#Interdependence - Western v Russia/China -----
model1 = model2 = model3 = model4 =NULL 
vcov_matrix1 = vcov_matrix2 = vcov_matrix3 = vcov_matrix4 = NULL 

formula1 <- shamed_dummy ~ 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ russ_china_shame_prior*usa + prior_no_russia_china

formula2 <- shamed_dummy ~ 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ russ_china_shame_prior*uk + prior_no_russia_china 

formula3 <- shamed_dummy ~ 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + russ_china_shame_prior*france + prior_no_russia_china 

formula4 <- shamed_dummy ~ 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + russ_china_shame_prior*western + prior_no_russia_china 

# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3, formula4)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))

data_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))
model4 <- glm(formula4, data = data_clean_all, family = binomial("logit"))


vcov_matrix1 <- vcovFWB(model1, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix4 <- vcovFWB(model4, cluster = data_clean_all$unsc_gwn, R = universal_r)


# Generate the table correctly
msummary(
  list(model1, model2, model3, model4), # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2,vcov_matrix3, vcov_matrix4), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c("*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  # Customize the bottom section of the table
  
  title= "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming}: Do Western Powers Respond to Russia and China?",
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/interdependence4.tex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')
)



#Controls for terrorism --------------------
formula1 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + t_a_lag
formula2 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + t_b_lag
formula3 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + k_a_lag
formula4 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + k_b_lag


# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3, formula4)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))

data_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

# This code should now run without error
model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))
model4 <- glm(formula4, data = data_clean_all, family = binomial("logit"))

vcov_matrix1 <- vcovFWB(model1, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix4 <- vcovFWB(model4, cluster = data_clean_all$unsc_gwn, R = universal_r)

msummary(
  list(model1, model2, model3, model4), # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2,vcov_matrix3, vcov_matrix4), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c("+" = .1, "*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  # Customize the bottom section of the table
  title =  "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming}: Includes Terrorism", 
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/appendix_terrorism.tex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')
)

#Controls for peace negotiations ------
formula1 <- shamed_dummy ~ ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior
formula2 <- shamed_dummy ~ best_fatality_estimate_lag_log*p5 +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior
formula3 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + neg_lag



# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))

data_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

# This code should now run without error
model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))

vcov_matrix1 <- vcovFWB(model1, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = data_clean_all$unsc_gwn, R = universal_r)

msummary(
  list(model1, model2, model3), # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2,vcov_matrix3), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c("+" = .1,"*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  # Customize the bottom section of the table
  title =  "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming}: Includes Peace Negotiations", 
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/appendix_negotiations.tex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')
)

#Includes UNSC member state support to NSAGs ------
formula1 <- shamed_dummy ~ ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior
formula2 <- shamed_dummy ~ best_fatality_estimate_lag_log*p5 +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior
formula3 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + ext_nsa_confirmed_lag

# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))

data_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

# This code should now run without error
model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))

vcov_matrix1 <- vcovFWB(model1, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = data_clean_all$unsc_gwn, R = universal_r)

# Generate the table correctly
msummary(
  list(model1, model2, model3), # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2,vcov_matrix3), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c("*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  # Customize the bottom section of the table
  title = "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming}: Include External Support to NSAGs", 
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/includes_ext_sup_to_nsag.tex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')
)




#Includes linear time trend ------
formula1 <- shamed_dummy ~ ext_confirmed_lag*p5  +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior
formula2 <- shamed_dummy ~ best_fatality_estimate_lag_log*p5  +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior
formula3 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag + year


# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))

data_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

# This code should now run without error
model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))

vcov_matrix1 <- vcovFWB(model1, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = data_clean_all$unsc_gwn, R = universal_r)

msummary(
  list(model1, model2, model3), # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2,vcov_matrix3), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c("*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  # Customize the bottom section of the table
  title =  "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming}: Includes Linear Time Trend", 
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/appendix_linear_time_trend.tex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')
)



####### Just P5 that look like E10  ------------------
formula0 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log+p5 +
  ext_confirmed_lag + battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag
formula1 <- shamed_dummy ~ ext_confirmed_lag*p5 
formula2 <- shamed_dummy ~ best_fatality_estimate_lag_log*p5 
formula3 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log+ext_confirmed_lag+ p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag
formula4 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*ext_confirmed_lag + p5 + battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag
formula5 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag


# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3, formula4)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))

data_clean_all <- shame2 %>%
  filter(!(.data$p5 == 1 & .data$total_govts_supporting > 7)) %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

# This code should now run without error
model0 <- glm(formula0, data = data_clean_all, family = binomial("logit"))
model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))
model4 <- glm(formula4, data = data_clean_all, family = binomial("logit"))
model5 <- glm(formula5, data = data_clean_all, family = binomial("logit"))

vcov_matrix0 <- vcovFWB(model0, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix1 <- vcovFWB(model1, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix4 <- vcovFWB(model4, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix5 <- vcovFWB(model5, cluster = data_clean_all$unsc_gwn, R = universal_r)

# Make sure you have a coefficient map for clean names
cm <- c(
  'ext_confirmed_lag' = 'External Support',
  'ext_confirmed_lag:p5' = 'External Support*P5',
  'p5:ext_confirmed_lag' = 'External Support*P5',
  'best_fatality_estimate_lag_log' = 'One-Sided Violence',
  'best_fatality_estimate_lag_log:p5' = 'One-Sided Violence*P5',
  'p5' = 'P5',
  'v2x_clphy_lag' = 'Physical Integrity',
  'battledeaths_lag_log' = 'Intensity',
  'incompatibility' = 'Secessionist',
  'pk_lag_agg_log' = 'Peacekeeping',
  'any_member_shamed_year_prior' = 'Prior Shame',
  'total_resolutions_lag_dummy' = 'Resolutions',
  'total_resolutions_log_lag' = 'Resolutions',
  '(Intercept)' = '(Intercept)')


# Generate the table correctly
msummary(
  list(model0, model1, model2, model3, model4, model5), # Use your fitted model object here
  vcov = list(vcov_matrix0, vcov_matrix1, vcov_matrix2,vcov_matrix3, vcov_matrix4, vcov_matrix5), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  #coef_map = cm,
  stars = c("+" = .1, "*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  # Customize the bottom section of the table
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/just_low_ext_sup_p5_A.tex",
  # Add table notes correctly
  notes = list('* p < 0.05, ** p < 0.01, *** p < 0.001',
               'Standard errors clustered by country are in parentheses.')
)

data_clean_all <- shame2 %>%
  filter(!(.data$p5 == 1 & .data$total_govts_supporting > 7)) %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  filter(year <2015) %>%
  na.omit()

# This code should now run without error
model0 <- glm(formula0, data = data_clean_all, family = binomial("logit"))
model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))
model4 <- glm(formula4, data = data_clean_all, family = binomial("logit"))
model5 <- glm(formula5, data = data_clean_all, family = binomial("logit"))

vcov_matrix0 <- vcovFWB(model0, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix1 <- vcovFWB(model1, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix4 <- vcovFWB(model4, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix5 <- vcovFWB(model5, cluster = data_clean_all$unsc_gwn, R = universal_r)

# Generate the table correctly
msummary(
  list(model0, model1, model2, model3, model4, model5), # Use your fitted model object here
  vcov = list(vcov_matrix0, vcov_matrix1, vcov_matrix2,vcov_matrix3, vcov_matrix4, vcov_matrix5), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  #coef_map = cm,
  stars = c("+" = .1, "*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/just_low_ext_sup_p5_B.tex",
  # Add table notes correctly
  notes = list('* p < 0.05, ** p < 0.01, *** p < 0.001',
               'Standard errors clustered by country are in parentheses.')
)


####### Split sample ------
formula1 <- shamed_dummy ~ any_member_shamed_year_prior + 
  ext_confirmed_lag + battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag
formula2 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log + battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag
formula3 <- shamed_dummy ~ any_member_shamed_year_prior + 
  best_fatality_estimate_lag_log+ext_confirmed_lag + battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag

# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn", "p5"))

e10_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  filter(p5==0)%>%
  na.omit()

p5_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  filter(p5==1)%>%
  na.omit()

# This code should now run without error
model1 <- glm(formula1, data = e10_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = e10_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = e10_clean_all, family = binomial("logit"))
model4 <- glm(formula1, data = p5_clean_all, family = binomial("logit"))
model5 <- glm(formula2, data = p5_clean_all, family = binomial("logit"))
model6 <- glm(formula3, data = p5_clean_all, family = binomial("logit"))

vcov_matrix1 <- vcovFWB(model1, cluster = e10_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = e10_clean_all$unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = e10_clean_all$unsc_gwn, R = universal_r)
vcov_matrix4 <- vcovFWB(model4, cluster = p5_clean_all$unsc_gwn, R = universal_r)
vcov_matrix5 <- vcovFWB(model5, cluster = p5_clean_all$unsc_gwn, R = universal_r)
vcov_matrix6 <- vcovFWB(model6, cluster = p5_clean_all$unsc_gwn, R = universal_r)


msummary(
  # The nested list for models is CORRECT for creating headers
  list(model1, model2, model3, model4, model5, model6), 
  #models_with_headers, # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2,vcov_matrix3, vcov_matrix4, vcov_matrix5, vcov_matrix6),
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c("+"=.1, "*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE,  
  title =  "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming}: Split Sample", 
  output = "latex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')) %>%
  add_header_above(header = c(" " = 1, "E10" = 3, "P5" = 3)) %>%
  save_kable(file = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/split_sample.tex")

###Intiator v. independent ------
formula1 <- shamed_dummy ~ e10_shamed_same_year_prior*p5 + p5_shamed_same_year_prior*p5 + 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag

formula2 <- shamed_dummy ~ e10_shamed_same_year_prior*best_fatality_estimate_lag_log*p5 + p5_shamed_same_year_prior*best_fatality_estimate_lag_log*p5 + 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag

formula3 <- shamed_dummy ~ e10_shamed_same_year_prior*ext_confirmed_lag*p5 + p5_shamed_same_year_prior*ext_confirmed_lag*p5 + 
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag+
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag

# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))


data_clean_all <- shame2 %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

# This code should now run without error
model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))


vcov_matrix1 <- vcovFWB(model1, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = data_clean_all$unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = data_clean_all$unsc_gwn, R = universal_r)

# Generate the table correctly
msummary(
  list(model1, model2, model3), # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2, vcov_matrix3), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c("*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  # Customize the bottom section of the table
  title = "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming} \\label{initiate2.tex}", 
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/initiate2.tex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')
)


# --- 1. Calculate Predicted Probabilities ---
# Create a grid of all combinations of the variables of interest.
# All other variables in model3 are held at their mean or mode.
pred_df <- predictions(
  model2,
  newdata = datagrid(
    p5 = c(0, 1), # 0 = E10, 1 = P5
    p5_shamed_same_year_prior = c(0, 1), # 0 = Not Shamed, 1 = Shamed
    best_fatality_estimate_lag_log = seq(
      min(shame2$best_fatality_estimate_lag_log, na.rm = TRUE),
      max(shame2$best_fatality_estimate_lag_log, na.rm = TRUE),
      length.out = 100
    )
  )
)

# --- 2. Create the Plot ---
# We map color to P5 status and line type to prior shaming status to get 4 lines.
plotz= ggplot(pred_df, aes(
  x = best_fatality_estimate_lag_log, 
  y = estimate, 
  color = factor(p5),
  linetype = factor(p5_shamed_same_year_prior)
)) +
  
  geom_line(linewidth = 1.1) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = factor(p5)), alpha = 0.15, linetype = "dotted") +
  
  labs(
    title = "",
    x = "One-Sided Violence (Logged Civilian Fatalities)",
    y = "Predicted Probability of Naming and Shaming",
    color = "Current Member Status",
    linetype = "Prior Action",
    fill = "Current Member Status"
  ) +
  
  scale_color_manual(
    values = c("0" = "#D55E00", "1" = "#0072B2"),
    labels = c("0" = "E10 Member", "1" = "P5 Member")
  ) +
  scale_fill_manual(
    values = c("0" = "#D55E00", "1" = "#0072B2"),
    labels = c("0" = "E10 Member", "1" = "P5 Member")
  ) +
  scale_linetype_manual(
    values = c("0" = "dashed", "1" = "solid"),
    labels = c("0" = "Not Previously Shamed by P5", "1" = "Previously Shamed by P5")
  ) +
  
  theme_bw() +
  theme(
    legend.position = "bottom", 
    legend.box = "vertical",
    # --- ADD THIS LINE TO MAKE THE LEGEND LINES LONGER ---
    legend.key.width = unit(1.5, "cm") 
  )
# 6. Save the plot to a JPEG file
ggsave(
  "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/sequence1.jpeg",
  plot = plotz,
  width = 8,
  height = 6,
  dpi = 300,
  bg = "white"
)


pred_df_e10 <- predictions(
  model2,
  newdata = datagrid(
    p5 = c(0, 1), 
    e10_shamed_same_year_prior = c(0, 1), # MODIFIED: Changed from p5_shamed...
    best_fatality_estimate_lag_log = seq(
      min(shame2$best_fatality_estimate_lag_log, na.rm = TRUE),
      max(shame2$best_fatality_estimate_lag_log, na.rm = TRUE),
      length.out = 100
    )
  )
)

# --- 2. Create the Plot ---
# The linetype aesthetic is now mapped to the E10 prior shaming variable
plotz_e10 = ggplot(pred_df_e10, aes(
  x = best_fatality_estimate_lag_log, 
  y = estimate, 
  color = factor(p5),
  linetype = factor(e10_shamed_same_year_prior) # MODIFIED: Changed from p5_shamed...
)) +
  
  geom_line(linewidth = 1.1) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = factor(p5)), alpha = 0.15, linetype = "dotted") +
  
  labs(
    title = "",
    x = "One-Sided Violence (Logged Civilian Fatalities)",
    y = "Predicted Probability of Naming and Shaming",
    color = "Current Member Status",
    linetype = "Prior Action",
    fill = "Current Member Status"
  ) +
  
  scale_color_manual(
    values = c("0" = "#D55E00", "1" = "#0072B2"),
    labels = c("0" = "E10 Member", "1" = "P5 Member")
  ) +
  scale_fill_manual(
    values = c("0" = "#D55E00", "1" = "#0072B2"),
    labels = c("0" = "E10 Member", "1" = "P5 Member")
  ) +
  # The legend labels are updated to reflect the new variable
  scale_linetype_manual(
    values = c("0" = "dashed", "1" = "solid"),
    labels = c("0" = "Not Previously Shamed by E10", "1" = "Previously Shamed by E10") # MODIFIED
  ) +
  
  theme_bw() +
  theme(
    legend.position = "bottom", 
    legend.box = "vertical",
    legend.key.width = unit(1.5, "cm") 
  )

# 6. Save the plot to a JPEG file
ggsave(
  "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/sequence2.jpeg", # MODIFIED: New filename
  plot = plotz_e10,
  width = 8,
  height = 6,
  dpi = 300,
  bg = "white"
)


pred_df_ext <- predictions(
  model3, # Make sure model3 contains the ext_confirmed_lag interaction
  newdata = datagrid(
    p5 = c(0, 1), 
    p5_shamed_same_year_prior = c(0, 1),
    ext_confirmed_lag = c(0, 1) # MODIFIED: Replaced the continuous OSV variable
  )
)

# --- 2. Create the Plot ---
# We now use geom_point and geom_errorbar for the categorical x-axis.
plot_ext = ggplot(pred_df_ext, aes(
  x = ext_confirmed_lag, 
  y = estimate, 
  color = factor(p5),
  linetype = factor(p5_shamed_same_year_prior)
)) +
  
  # MODIFIED: Changed from geom_line to geom_point and geom_errorbar
  geom_point(position = position_dodge(width = 0.4), size = 3) +
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high),
    width = 0.2, 
    position = position_dodge(width = 0.4)
  ) +
  
  # MODIFIED: Updated labels for the new variable
  labs(
    title = "",
    x = "External Support (Lagged)",
    y = "Predicted Probability of Naming and Shaming",
    color = "Member Status",
    linetype = "Prior Action"
  ) +
  
  # MODIFIED: Updated the x-axis scale for the new dummy variable
  scale_x_continuous(
    breaks = c(0, 1),
    labels = c("No External Support", "External Support")
  ) +
  
  scale_color_manual(
    values = c("0" = "#D55E00", "1" = "#0072B2"),
    labels = c("0" = "E10 Member", "1" = "P5 Member")
  ) +
  scale_linetype_manual(
    values = c("0" = "dashed", "1" = "solid"),
    labels = c("0" = "Not Previously Shamed by P5", "1" = "Previously Shamed by P5")
  ) +
  
  theme_bw() +
  theme(
    legend.position = "bottom", 
    legend.box = "vertical",
    legend.key.width = unit(1.5, "cm") 
  )

# 6. Save the plot to a JPEG file
ggsave(
  "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/sequence3.jpeg", # MODIFIED: New filename
  plot = plot_ext,
  width = 8,
  height = 6,
  dpi = 300,
  bg = "white"
)

pred_df_e10 <- predictions(
  model3, # Make sure model3 contains the e10_shamed_same_year_prior interaction
  newdata = datagrid(
    p5 = c(0, 1), 
    e10_shamed_same_year_prior = c(0, 1), # MODIFIED: Replaced the p5_shamed... variable
    ext_confirmed_lag = c(0, 1)
  )
)

# --- 2. Create the Plot ---
# The linetype is now mapped to 'e10_shamed_same_year_prior'.
plot_e10 = ggplot(pred_df_e10, aes(
  x = ext_confirmed_lag, 
  y = estimate, 
  color = factor(p5),
  linetype = factor(e10_shamed_same_year_prior) # MODIFIED: Replaced the p5_shamed... variable
)) +
  
  geom_point(position = position_dodge(width = 0.4), size = 3) +
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high),
    width = 0.2, 
    position = position_dodge(width = 0.4)
  ) +
  
  # MODIFIED: Updated labels for the new variable
  labs(
    title = "",
    x = "External Support (Lagged)",
    y = "Predicted Probability of Naming and Shaming",
    color = "Member Status",
    linetype = "Prior Action"
  ) +
  
  scale_x_continuous(
    breaks = c(0, 1),
    labels = c("No External Support", "External Support")
  ) +
  
  scale_color_manual(
    values = c("0" = "#D55E00", "1" = "#0072B2"),
    labels = c("0" = "E10 Member", "1" = "P5 Member")
  ) +
  
  # MODIFIED: Updated the legend labels
  scale_linetype_manual(
    values = c("0" = "dashed", "1" = "solid"),
    labels = c("0" = "Not Previously Shamed by E10", "1" = "Previously Shamed by E10")
  ) +
  
  theme_bw() +
  theme(
    legend.position = "bottom", 
    legend.box = "vertical",
    legend.key.width = unit(1.5, "cm") 
  )

# 6. Save the plot to a JPEG file
ggsave(
  "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/sequence4.jpeg", # MODIFIED: New filename
  plot = plot_e10,
  width = 8,
  height = 6,
  dpi = 300,
  bg = "white"
)

######Matched Sample ------
external = read.csv("/Users/cbrandt/Dropbox/Data/UCDP/ucdp-esd-ty-181.csv")
external = external[which(external$actor_nonstate==0&external$ext_nonstate==0),]
external2 = external %>% 
  group_by(ext_name, actor_id) %>%
  arrange(year) %>%
  mutate(
    gwno_a = countrycode(ext_name, "country.name", "gwn"),
    ext_confirmed_lag = lag(ext_sup, n = 1)
  )  %>%
  ungroup()%>%
  group_by(gwno_a, year)%>%
  summarise(total_govts_supporting=sum(ext_confirmed_lag))%>%
  ungroup() 


security_council <- read.csv("/Users/cbrandt/Dropbox/AAA_Projects/UNSC/data/UNSCdata.csv")
security_council = security_council %>% 
  mutate(
    year = as.numeric(year),
    p5 =  factor(
      p5,
      levels = c("1", "0"),
      labels = c("P5", "E10")
    ),
    gwno_a = if_else(
      country == "St. Vincent",
      57,
      countrycode(country, "country.name", "gwn")
    )
  ) %>% 
  left_join(external2, by=c("year", "gwno_a")) %>%
  mutate(total_govts_supporting = ifelse(is.na(total_govts_supporting)==TRUE, 0, total_govts_supporting))


# --- Step 1: Bin the Continuous Data into Categories (Same as before) ---
plot_data <- security_council %>%
  mutate(
    support_category = case_when(
      total_govts_supporting == 0 ~ "0",
      total_govts_supporting == 1 ~ "1",
      total_govts_supporting == 2 ~ "2",
      total_govts_supporting >= 3 & total_govts_supporting <= 5 ~ "3-5",
      total_govts_supporting > 5 ~ "6+",
      TRUE ~ "Other"
    )
  ) %>%
  mutate(
    support_category = factor(support_category, levels = c("0", "1", "2", "3-5", "6+"))
  ) %>%
  count(p5, support_category, name = "frequency")

# --- Step 2: Create a Single, Dodged Bar Chart with Custom Colors ---
final_plot_single <- ggplot(plot_data, aes(x = support_category, y = frequency, fill = p5)) +
  # Use position="dodge" to place bars side-by-side
  geom_col(position = "dodge") +
  
  # Add custom colors for P5 and E10
  scale_fill_manual(values = c("E10" = "#D55E00", "P5" = "#0072B2")) +
  
  # Add the axis break (adjust values as needed for the combined data)
  scale_y_break(breaks = c(80, 180)) +
  
  # Add labels and styling
  labs(
    title = "",
    x = "Number of Governments UNSC Member States are Providing Support to, Yearly",
    y = "Frequency",
    fill = "Council Bloc" # Renames the legend title
  ) +
  theme_bw() +
  theme(
    plot.title = element_text(hjust = 0.5),
    legend.position = "bottom"
  )

# --- Save the Final Plot ---
ggsave(
  "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/plots/number_of_states_providing_aid_to_governments.jpg", 
  plot = final_plot_single, 
  width = 10, 
  height = 7, 
  dpi = 300
)
# --- (Your code to write the CSV remains the same) ---
shame2$total_govts_supporting = log(shame2$total_govts_supporting+1 )

match_no_replace <- matchit(
  p5 ~ total_govts_supporting, 
  data = shame2,
  method = "nearest", 
  replace = FALSE,
  caliper = .1
  # This line drops treated units that don't have a close match
)

summary(match_no_replace)

matched_data_no_replace <- match.data(match_no_replace)

t.test(total_govts_supporting ~ p5, data = matched_data_no_replace)


#######Main Analysis ======
formula1 <- shamed_dummy ~ ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag +
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag+ any_member_shamed_year_prior

formula2 <- shamed_dummy ~ any_member_shamed_year_prior +
  best_fatality_estimate_lag_log*p5  +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag

formula3 <- shamed_dummy ~ any_member_shamed_year_prior +
  best_fatality_estimate_lag_log*p5 + ext_confirmed_lag*p5 +battledeaths_lag_log + v2x_clphy_lag + 
  incompatibility + pk_lag_agg_log + total_resolutions_log_lag

# 2. Put all formulas into a list for easy processing
formula_list <- list(formula1, formula2, formula3)

# 3. Dynamically extract all unique variable names from the formulas
vars_from_formulas <- unique(unlist(lapply(formula_list, all.vars)))

# 5. Create the final, complete list of variables
all_needed_vars <- unique(c(vars_from_formulas, "year", "unsc_gwn"))

data_clean_all <- matched_data_no_replace %>%
  dplyr::select(all_of(all_needed_vars)) %>%
  na.omit()

# This code should now run without error
model1 <- glm(formula1, data = data_clean_all, family = binomial("logit"))
model2 <- glm(formula2, data = data_clean_all, family = binomial("logit"))
model3 <- glm(formula3, data = data_clean_all, family = binomial("logit"))

vcov_matrix1 <- vcovFWB(model1, cluster = ~unsc_gwn, R = universal_r)
vcov_matrix2 <- vcovFWB(model2, cluster = ~unsc_gwn, R = universal_r)
vcov_matrix3 <- vcovFWB(model3, cluster = ~unsc_gwn, R = universal_r)

# Generate the table correctly
msummary(
  list(model1, model2, model3), # Use your fitted model object here
  vcov = list(vcov_matrix1, vcov_matrix2,vcov_matrix3), # Use your calculated vcov matrix
  
  # These arguments produce a clean, publication-quality table
  booktabs = TRUE,
  coef_map = cm,
  stars = c("+" = .1, "*" = .05, "**" = .01, "***" = 0.001),
  exponentiate = TRUE, 
  # Customize the bottom section of the table
  title = "Odds Ratios from Logistic Regression on \\textit{Naming and Shaming}  \\label{match_no_replacement}", 
  output = "/Users/cbrandt/Dropbox/Apps/Overleaf/Naming and Shaming Rebel Groups/regression_analysis2/matched_sample.tex",
  # Add table notes correctly
  notes = list('Bootstrapped standard errors clustered by country.')
)

